arXiv: 1504.04594vl [q-fin.MF] 17 Apr 2015 


A Posteriori Error Estimator for a Front-Fixing 
Finite Difference Scheme for American Options 

Riccardo Fazio 


Abstract —For the numerical solution of the American option 
valuation problem, we provide a script written in MATLAB 
implementing an explicit finite difference scheme. Our main 
contribute is the definition of a posteriori error estimator for 
the American options pricing which is based on Richardson’s 
extrapolation theory. This error estimator allows us to find a 
suitable grid where the computed solution, both the option price 
field variable and the free boundary position, verify a prefixed 
error tolerance. 

Index Terms —American options, free boundary problem, fi¬ 
nite difference scheme, Richardson’s extrapolation, a posteriori 
error estimator. 


I. Introduction 

In the market for financial derivatives, the most important 
problem is the so-called option valuation problem or in a 
few words: the problem of computing a fair price for a 
given option. An American call (put) option is a contract 
written on an underlying asset and gives the holder the right 
to buy (sell) the asset for a pre-specified price or strike price 
on or before a pre-specified date also indicated as maturity. 
Unlike European options where the holder can exercise the 
option only on the maturity date, the possibility of early 
exercise makes the pricing of American options a problem 
in stochastic optimization. While a closed-form solution for 
the price of European options is derived in the celebrated 
work by Black and Scholes |1] and by Merton [5], there 
exists no analogous result for American options. The reason 
can be explained as follows: while the governing differential 
equation is still the one obtained by Black and Scholes 
(T), McKean 0 and Merton (4) show that the price of 
an American option satisfies boundary conditions governed 
by a boundary that is not known a priori and needs to be 
computed as part of the solution itself. Such problems are 
called free-boundary problems. In particular, the American 
call option problem is a free boundary problem defined on a 
finite interval. On the other hand, the American put option 
problem is a free boundary problem defined on a semi- 
infinite interval so that it is a non-linear problem complicated 
by a boundary condition at infinity. Therefore, both such 
derivatives of financial markets must be priced by numerical 
or analytical approximations. 

Within analytical approximations, MacMillan (5| and 
Barone-Adesi and Whaley (61 defined a quadratic approx¬ 
imation for American put option. These methods are not 
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convergent and have trouble pricing long-maturity options 
accurately. To correct this problem, Ju and Zhong (1999) 
develop an approximation based on the method proposed 
in Barone-Adesi and Whaley. While this improved method 
prices long-maturity options more accurately, it is still not 
convergent. Johnson (8| used an interpolation scheme to 
price the American put option, Geske and Johnson [9j 
derived a valuation formula expressed in terms of a series of 
compound-option functions for the same reason, Bunch and 
Johnson m propose a modified two-point Geske-Johnson 
approach. Carr and Faguet GD view put options as the limit 
of a sequence of perpetual option values which are subject 
to default risk, and use this view to deriving approximations 
for the price of an option. More recently, Zhu (12) derives 
a semi-closed form solution for the price of the option as 
a Taylor series expansion consisting of infinite terms, but 
requiring thirty terms for an accurate option price. 

Within numerical approximations, the most popular meth¬ 
ods for pricing American option can be classified to lat¬ 
tice method, Monte Carlo simulation, and Finite difference 
method. The lattice method was first introduced by Cox 
et al. 03) , and its convergence was proved by Amin and 
Khanna~]l4|. Fu (B), og applied Monte Carlo method 
along with gradient-based optimization techniques to price 
American style options. The application of a finite difference 
method to price American options was initiated by Brennan 
and Schwartz 07], 00 and Schwartz GD- Jaillet et al. (D 
proved the finite difference method convergence. A front- 
fixing finite difference method was proposed by Wu and 
Kwok (1997), Nielsen et al. pi] , and Company et al. (22) 
to compute option prices. The front-fixing method utilizes a 
change of variables to transform the free boundary problem 
into a nonlinear problem on a fixed domain. Nielsen et al. 
GD also propose a penalty method to price American put 
options, where the unknown boundary is removed by adding 
a penalty term, again leading to a nonlinear problem posed 
on a fixed domain. 

In this paper, we list a script written in MATFAB imple¬ 
menting a finite difference scheme for the numerical solution 
of the American option models of financial markets. We 
implemented both the method defined by Company et al. 
[22 1 and the method developed by Nielsen et al. GD and 
found that the first method implementation is more efficient 
than the second one. Our main contribute is the definition of 
a posteriori error estimator for the American options pricing 
which is based on Richardson’s extrapolation theory. This 
error estimator allows us to find a suitable grid where the 
computed solution, both the option price field variable and 
the free boundary position, verify a prefixed error tolerance. 


II. American put option 


Let us suppose that at time t the price of a given underlying 
asset is S. We consider here the following mathematical 
model for the value P = P(S, r) of an American put option 
to sell the asset: 


9P = 1 2.2 &P 

dr 2° dS 2 


dP 

dS 


P(S, 0) = ma x(E — S, 0) 


lim P(S , r) = 0 , 

S—> oo 


rP , 

S*{0) = E , 


p(s*(t),t) = e-s*(t) , og{sr(T),T) 
P(S, t) = E — S , 0 <S<S*(t) 


(1) 


- 1 , 


where r = T—t, t denotes the time to maturity T, S*(t) is a 
free boundary, that is the unknown early exercise boundary, 
cr, r and E are given constant parameters representing the 
volatility of the underlying asset, the interest rate and the 
exercise price of the option, respectively, and the governing 
equation is defined on 0 < r < T, S*(t) < S < oo. 

To fix the free boundary we apply the dimensionless new 
variables 


x = In 


Sf = 

p(x,r) = 


s f( r ) ’ 

S*(r) 

E ’ 

P(xSf{T),T) 
E 


( 2 ) 


see Wu and Kwok ]23j. From the variables transformation 
defined by <[2j) follows that Sj(t) is mapped on the fixed line 
x = 0, 0 < p{x,r) < 1 and 0 < S/(t) < 1. By using 0, 
the put option problem 0 can be rewritten as follows 


dp = 1 2 d2 P 

dr 2° dx 2 


r — 


p(x, 0) = 0 for 0 < x , 
lim p(x, t) = 0 , 

x—too 

p(0,r) = 1 - S f (r) , 


dp 

77 - rp 

ox 


Sf( 0 ) = 1 , 


dp 

dx 


(0,7-) = S f ( 7 -) 


1 dSf . .dp 
Sf(r) dr dx 1 

(3) 

(4) 

(5) 

( 6 ) 


that has to be solved on the domain defined by 0 < t < T 

and 0 < x < 00 . 


III. An explicit finite difference scheme 

To solve the problem ([3]j6]» numerically, we introduce a 
truncated boundary x^. which is a suitable large value where 
it would be convenient to impose the asymptotic boundary 
condition. In other words, we replace the asymptotic bound¬ 
ary condition 0 with the side condition 

p{xoo,t) = 0. (7) 

For the choice of Xoo and the accuracy of the related 
numerical solution, we can refer to the study by Kangro and 
Nicolaides f24 ). On the other hand, the boundary condition 
at infinity can be enforced exactly by using a non-standard 
finite difference scheme, and this has been shown for the 
numerical solution of the so-called perpetual American put 
option in [25]. 


Next, by setting an integer J and a positive value //, we 
define the step-sizes 


the integer N 


*£00 

At = pAx 2 , 

(8) 

N = 

1 

(9) 


where [■] : 1R + —> N is the ceiling function which maps a 
real number to the least integer that is greater than or equal 
to that number. Therefore, /< is the grid-ratio 


At 

^ Ax 2 


( 10 ) 


So that, within the finite domain, we can introduce a mesh 
of grid-points 


Xj = j Ax , t n = nAt , (11) 

for j = 0,1,..., J and n = 0,1,..., N. 

We would like to define a numerical scheme that allows 
us to compute the grid values 


Pj ~ P(xj,t n ) , (12) 

for j = 0,1,..., J and n = 0,1,..., IV — 1, N and the free 
boundary values 

S] « Sf(t n ) , (13) 


for n = 0,1,..., N — 1, N. To this end let us consider the 
explicit finite difference scheme 


Po 


n+1 


= 1 2 pU- 2 p]+p7+i 

At 2 (A.t) 2 

P? + i-Pj-i 


+ r —— 


2Ax 


+ 


(14) 


+ 


cin+t on n 

b f Pj+l 


Pj-1 


A tS n f 


2Ax 


r Pj 


for j = 1, 2,..., J — 1 and n = 0,1,..., N — 1. For our 
specific problem p" and S'f are given and our goal is to 
compute p" +1 and ,5')' + l . If we apply some simple algebra, 
then we can rewrite the explicit finite difference scheme as 


ora+l on n _ n 

„«+! _ h n n rn n _ L _/ Pj+1 Pj-1 

Pj a Pj -1 + °Pj + c Pj +1 + AtS n 2 Ax 

(15) 

for j = 2,3,..., J — 1 and n = 0,1,..., N — 1, where 

. 2 ' 


P 

a= 2 


a 1 — r —— Ax 


6 = 1 — pa 2 — i'At 


( 16 ) 


C 2 


cr" — | r + — ) Ax 


Now, we have to take into account the side conditions. From 
the two initial conditions 0- we obtain 

P° = 0 , S°f = l, (17) 


for j = 0,1,..., J. From the boundary conditions (|7j, we 
get 


Pj = 0 , 


(18) 





















tor n = 0,1,.... A. From the two boundary conditions (6), 
using a central finite difference formula, we derive 


Pi - P-i 
2Ax 


= -S' 


m n _ i on 

/ > Po - 1 ~ J 


(19) 


where X-\ = —Ax is a fictitious point out of the com¬ 
putational domain. Moreover, by considering the governing 
differential equation (3) at Xo = 0, r > 0 and taking into 
account the side conditions (6) one gets a new boundary 
condition: 


cr“ d 2 p , a 


T§^ (oV) + t Si{t) ~ t = 0 


a 2 p n _ l -2p™+p™ a 


2 S ?-r = 0. 


2 Ax 2 

Now, we can eliminate the value of p" 1 from equations (19 1 


and (21 1 to get 


rAx 2 f Ax 2 \ 

Pi = 1 + —^-( 1 + Ax + ) S n f 


S™ +1 = d n S] 


for n = 0,1,..., A — 1, where 
rAx 2 


1 + 


d n = 


- [apo + bp™ + cp% - 


P 2 - P Q 
2Ax 


P'i - Po 
2Ax 


+ 1 


Ax 2 ' 

Ax+ — ) S ™ f 


compute S^ +1 according to (|23|l and apply the free 


5) for n = 0,1,... N — 1, compute d n according to (24 1 , 
compute S^ +1 acco 
boundary conditions 


IV. Positivity, monotonicity and stability 

In this section, we recall the theoretical results that make 
the explicit difference scheme suitable for doing numerical 
studies of the American put option. For the sake of simplicity, 
let us define the numerical solution vector at time level t n 

as p" = (pff, p^,...,p n j) T . 

Lemma 1: (Due to Company et al. (22).) If Af and Ax 
verify the two inequalities 


Ax < 


( 20 ) 


|r — rx 2 / 21 


Af , r / er 2 /2 


see Wu and Kwok (23) , Zhang and Zhu (26) or Kwok (27] 
p. 341] , and its central finite difference discretization 


At < 


Ax 2 


/t-2 


a 2 + rAx 2 


(25) 


(26) 


( 21 ) 


( 22 ) 


If we evaluate the numerical scheme <0 for j = 1 and take 
into account (221 for n = n+ 1, then the free boundary S™ +1 
can be definecThy 


(23) 


(24) 


We are now ready to define the algorithm: 

1) Input c, r, E, T, J, p and x^; 

2) Define the grid (Xj, Ax, t n , Af) according to equations 
(8) and (IT) ; 

3) for j = 0,1, 2,..., J do pj = 0 end, set S ° = 1; 

4) compute a, b and c according to (]T6|); 


then the coefficients a, b and c are non-negative. If r = cr 2 /2, 
then the non-negativity of these coefficients is verified under 
the condition (26). 

These two inequalities follow from positivity preserving 
considerations related to the explicit finite difference scheme. 
The theoretical framework is similar to the one used by 
Friedrichs (28) to study positivity preserving finite differ¬ 
ence schemes for the advection equation, see also Fazio 
and lannelli 29 . Moreover, the coefficients of p" + l , for 
i = j — 1 ,j,j + 1, in our difference scheme © might be 
regarded, in the limit for Af —> 0, as probability values. In 
fact, they sum up to 1 — r At and by imposing their positivity 
we get the two inequalities (25)-(|26). 

Theorem 1: (Due to Company et al. 122].) Let {p™, } 
be the computed numerical solution antlT/' be defined by 
equations (24) , then under hypothesis of the Lemma |T] for 
sufficiently small values of Ax, we have: 

• {S'j} is positive and non-increasing monotone for n = 

0,1,..., JV; 

• the vectors p n have positive components for n = 
0,1,...,A; 

• the vectors p n are non-increasing monotone with respect 
to the j for each fixed n = 0,1,..., N. 

As far as the stability of the explicit finite difference 
scheme is concerned we can introduce the definition: the 
numerical scheme is said to be || • Hoo-stable if, for every 
mesh in the computational domain [0, Xqo] x [0, X 1 ], there 
exists a positive constant C such that 


p™+ 1 = l - S n f +1 


P n i +1 = 1 + 


rAx 2 


- 1 + Ax + 


Ax 2 


cm+l 


compute 


a = a — 


(™+l _ gn 

2AxS] 


gn+1 _ Q*n 

2A xSf 


and for j = 2.1.— 1 compute the values of p" 4 
according to 


p] +1 =a n p]_i+bp]+c n p; +1 . 

end, set p " = 0, end. 

The implementation of this algorithm in MATLAB has been 
used to get the numerical results reported below. The script 
file is listed in the Appendix. 


lbloo<C for n = 0,1,..., A , (27) 

where C is independent on Af, Ax and n, see Company et 
al. (30). 

Theorem 2: (Due to Company et al. (22).) Under the 
hypothesis of Theorem[I]the explicit finite difference scheme 
for the fixed domain problem (|3]|6]) is || • Hoc-stable. 

V. A POSTERIORI ERROR ESTIMATOR 

For a scalar U of interest, either a value of the solution 
p™ or a free boundary value S'j , the numerical error e can 
be defined by 

e = u — U , (28) 

where u is the exact, usually unknown, value. When the 
numerical error is caused prevalently by the discretization 
error and in the case of smooth enough solutions the global 
























error can be decomposed into a sum of powers of the inverse 
of N 

PO / ^ \ Pi / 1 \ P2 


u — Un + Cq ( jy 


N 


Ci - +C 2 - 


1 


N 


+ ■ 


(29) 

where Co, Ci, C 2 , ... are coefficient that depend on u and 
its derivatives, but are independent on N, and po, p \, pi, ■ .. 
are the true orders of the discretization error, see Schneider 
and Marchi (31) and the references quoted therein. Each pt 
is usually a positive integer with po < pi < P 2 < ■ ■ ■ and 
all together they constitute an arithmetic progression of ratio 
Pi —po- The value of po is called the asymptotic order or the 
order of accuracy of the method or of the numerical value 
U N . By replacing into equation ( |29| N = N g and N = N g+ 1 
and subtracting, to the second obtained equation the first 
times (l/q) Po , q = N g+:i /N g , we get the first extrapolation 
formula 


n , £4+i - u g 

U 9 +1 + 


(30) 


qPo - 1 

that has a leading order of accuracy equal to p\. This type 
of extrapolation is due to Richardson (32), (33). Taking 


into account equation (30i we can conclude that the error 
estimator by a first Richardson’s extrapolation is given by 


£4+i - £4 


(31) 


qP “-I 

where po is the order of the numerical method used to 
compute the numerical solutions. Hence, (3l) gives the real 
value of the numerical solution error without knowledge of 
the exact solution. In comparison with © a safer error 
estimator can be defined by 


e s — £4+i ££c 


(32) 


Po 


(33) 


Of course, po can be estimated with the formula 

log(|£/ g ~u|)-log(|£/ g +i -m|) 
log(g) 

where u is again the exact solution (or, if the exact solution 
is unknown, a reference solution computed with a suitable 
large value of N), and both u and U g +'i are evaluated at the 
same grid-points of U g . 

Within the above framework, in order to improve the 
numerical accuracy by using only a small number of grid- 
nodes, we can generalize © introducing the following 
repeated extrapolation formula 


U, 


g+l,k+l 


= u, 


u, 


3+1, k 


9+1, k 


,-u, 


g,k 


(34) 


qPk - l 

where g £ {0,1, 2,..., G - 1}, k £ {0,1, 2,..., G - 1}, 
q = N g+ i/N g is the grid refinement ratio, and p^ is the 
true order of the discretization error. The formula © is 
asymptotically exact in the limit as Nq goes to infinity if 
we use uniform grids. We notice that to obtain each value of 
U g + 1 ,/c+i requires two computed solutions U in two adjacent 
grids, namely g + 1 and g at the extrapolation level /,:. For 
any g, the level k = 0 represents the numerical solution of 
U without any extrapolation. We recall that the theoretical 
orders of accuracy of the numerical values U gt k, with N = 
N g and k extrapolations, verify the relation 


Pfc = Po + k(pi -p 0 ) , 


(35) 


where this equation is valid for k £ {0,1, 2,..., G — 1}. 



X 


Fig. 1. Instability for N = 100, fi = 27 and Xqo = 1. 


VI. Numerical results 

We introduce a numerical test for the finite difference 
schemes defined below. To this end, we consider the Ameri¬ 
can put option problem (3]|6) with the following parameters: 

r = 0.1 , £7 = 0.2, E = T=1. (36) 

First of all, we are interested to define a suitable value of 
the truncated boundary x^. Then, we investigate numerically 
how the choice of the value of Xoo influences the numerical 
solution. To take a simple approach, we can monitor the final 
free boundary computed values Sj . Some sample results are 
reported in Table [I] and indicate that we can set = 1. 

TABLE I 

Free boundary location at t = T. 


Xoo 

N = 10 

N = 20 

O 

II 

i 

0.8676354534435 

0.8655750222427 

0.86438663362975 

2 

0.8676354534435 

0.8655750222427 

0.86438663362975 

4 

0.8676354534435 

0.8655750222427 

0.86438663362975 


Fig. □ shows an unstable computation. For the sake of 
clarity, in Fig. [I] we display only the initial condition and 
the numerical solution at t = T, that is the last time iterate. 
As usual, the instability manifest itself with large oscillations 
between positive and negative values. Therefore, as far as 
our explicit finite difference scheme is concerned, we have to 
find a compromise between accuracy and stability. Of course, 
from now on, the chosen grid-spacings are defined in order 
to verify the stability inequalities 


_-<(26 l 

From the results listed in Table [I] we realize that, fixed a 
value of the truncated boundary Xoo, the computed values of 
Sj 1 for different values of the grid-steps are in agreement 
only for the first two decimal places. Then, we decided 
to improve the numerical accuracy by performing a mesh 
refinement. Moreover, we applied repeated Richardson’s ex¬ 
trapolation to improve the numerical accuracy. Fet us recall 
that the explicit difference scheme is first order in time and 
second order in space both for the field variable and the 
free boundary value, i.e. the truncation error is of the order 
0(Ai 2 ) + 0(AtAx 2 ). We will use this result below when 





























we perform a mesh refinement keeping constant the grid- 
ratio, i.e. /i, so that we end up with second-order truncation 
error T" = 0(Af 2 ) in time and, therefore, the global error 
is first-order, that is the p (J value defined above is equal to 
one. We remark that in our case the sequence of q Pk , for 
k = 0,1,, is given by 4, 16, 64, 256, 1024, ..., that 

is q = 4 and p/~ = k + 1, for k = 0,1,_In Table [II] 

we report sample numerical results for the benchmark value 
Sj . Within the same Table, we report the results obtained 


TABLE II 

Richardson’s repeated extrapolations for the free boundary 

VALUE AT t = T. . 


N 

U g , 0 

U g , 1 

Ug, 2 

CO 

UgA 

5 

0.871621 





20 

0.865575 

0.863560 




80 

0.863700 

0.863075 

0.863043 



320 

0.863071 

0.862861 

0.862847 

0.862844 


1280 

0.862859 

0.862788 

0.862783 

0.862782 

0.862782 

5120 

0.862788 

0.862764 

0.862763 

0.862762 

0.862762 


by repeated Richardson’s extrapolations. We remark that the 
next extrapolated value from Table [II] is Upp = 0.862762 so 

0.862762. This value can 
^ ss 0.86222 computed by 
0.8628 found by Company et 
al. 


that our benchmark value is SjY ~ 
be compared with the values S^ 
Nielsen et al. [21] and Sj7 


Next, we indicate how to use the error estimator defined 
by equation ©, or alternatively by equation ( |32| ). Let us 
assume that our goal is to solve the American put option 
problem with a given tolerance e, where 0 < e < 1. To 
this end we should solve the given problem twice, for two 
grids defined with given values of J g = J and J g +i = 2J 
of space intervals but for the same value of the grid-ratio //, 
The corresponding time intervals N g and N g+ 1 verify the 
relation q = N g+ i/N g . Hence we can apply (component¬ 
wise) to p n and to SJ the error estimator formula (31 1 , or 
I. Then, we can verify whether, for n = 1,2,... . A, 


\\e r (p n )\\oo < e \e r (SJ)\<e. (37) 


If the two inequalities (37 1 hold true, for n = 1,2 ,N, 
then we can accept the numerical solution computed on 
the grid defined by J g +i and A r , ;+1 , otherwise we have to 
increase these two integers and repeat the calculations. 


Fig. 0 shows the error estimator results computed by 
setting e = 0.005. We fixed p = 20 and started with Jo = 5 
and Ji = 10 repeating the computation by doubling the 
number of spatial grid-intervals if the required accuracy was 
not achieved. Our algorithm stopped when J 4 = 80 that for 
p = 20 means N± = 320. 

For the sake of completeness, in Fig.|3]we plot pj r versus 
xj and SJ versus t n , these results were obtained by the 
finite difference scheme with N = 320 and p = 20. From 
the numerical results shown in figure [ 2 ] we can easily realize 
that the greatest errors are found within a few time steps. This 
suggests that a better accuracy can be achieved, without the 
use of Richardson’s extrapolation, by developing an adaptive 
version of the explicit finite difference scheme. 



Fig. 2. Numerical estimated errors e, r (p n ) und e r (SJ) versus t n . 



S f n of APO PDE with ji= 20 



Fig. 3. Numerical results: top frame pj versus Xj, bottom frame SJ 
versus t n - 


Appendix 

The MATLAB script file 

Here we list the basic algorithm written in MATLAB. The 
script file is called APOef ds . m and it can be easily modified 
to apply our a posteriori error estimator defined by equation 






















% APOefds.m 

% APO equation on [0,xinf] by forward 
% in time and central in space 
% J intervals in x, N intervals in t 
% mi = k/h"2, k time step, h space step 
clear all; help APOefds 
r = .1; % parameters 

T = 1; % by Nielsen et al. 2002 

sigma = .2; 

E - 1; 
xinf = 1; 

J - 80; 
mi = 2 0 
h - xinf/J 

k = mi*h"2 

N = ceil(T/k) 

x = 0:h:xinf; 

Stabl =sigma"2/(abs(r-.5*sigma"2))-h 
Stab2 =h"2/(sigma"2+r*h"2)- k 
if (StablcO | Stab2<0), break, end 
w = zeros(J+l,1); 
s(1) = 1; 
t(1)=0; 

A = .5*mi*(sigma"2-(r-.5*sigma"2)*h); 

B = 1-mi*sigma"2-r*k; 

C = .5*mi*(sigma"2+(r-.5*sigma"2)*h); 

Al = l+r*h"2/sigma"2; 

B1 = l+h+.5*h"2; 

for n=l:N %time loop 

dpi = .5* (w (3, 1) -w (1,1) ) /h; 

D = (Al-(A*w(1,1)+B*w(2,1)+ ... 

C*w(3,1)-dpi))/(dpl+Bl*s(n)); 
s (n+1) = D*s (n); 

if (s(n+l)<0 I s(n+l)>l), break, end 
w (1,2) = 1 — s (n+1); 
w(2,2) = Al-Bl*s(n+1); 

AM = A-.5*(s(n+1)-s(n))/(h*s(n)); 

CM = C+.5*(s (n+1)-s(n))/(h*s(n)); 
for j = 3:J %space loop 

w(j,2) = AM*w(j-l,l)+ ... 

B*w(j,1)+CM*w(j+l,1); 
end 

w(J+l,2) = 0; 
w ( : , 1) = w ( : , 2) ; 
t(n+1) = t(n)+k; 

end 

plot(t, s,'LineWidth',2.5) 
axis([0 T .86 1]); grid 
xlabel(' t' ) 

ylabel('S_f"n \approx S_f (t_n)') 
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